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PIECEWISE SMOOTH SYSTEMS NEAR A 
CO-DIMENSION 2 DISCONTINUITY MANIFOLD: 
CAN ONE SAY WHAT SHOULD HAPPEN? 

LUCA DIECI AND CINZIA ELIA 


Abstract. We consider a piecewise smooth system in the neighborhood of a co-dimension 
2 discontinuity manifold E (intersection of two co-dimension 1 manifolds). Within the 
class of Filippov solutions, if E is attractive, one should expect solution trajectories to 
slide on E. It is well known, however, that the classical Filippov convexification method¬ 
ology does not render a uniquely defined sliding vector field on E. The situation is further 
complicated by the possibility that, regardless of how sliding on E is taking place, during 
sliding motion a trajectory encounters so-called generic first order exit points, where E 
ceases to be attractive. 

In this work, we attempt to understand what behavior one should expect of a solution 
trajectory near E when E is attractive, what to expect when E ceases to be attractive 
(at least, at generic exit points), and finally we also contrast and compare the behavior 
of some regularizations proposed in the literature, whereby the original piecewise smooth 
system is replaced -in a neighborhood of E- by a smooth differential system. 

Through analysis and experiments in R® and we will confirm some known facts, 
and provide some important insight: (i) when E is attractive, a solution trajectory 
indeed does remain near E, viz. sliding on E is an appropriate idealization (of course, 
in general, one cannot predict which sliding vector field should be selected); (ii) when 
E loses attractivity (at first order exit conditions), a typical solution trajectory leaves 
a neighborhood of E; (iii) there is no obvious way to regularize the system so that the 
regularized trajectory will remain near E as long as E is attractive, and so that it will 
be leaving (a neighborhood of) E when E looses attractivity. 

We reach the above conclusions by considering exclusively the given piecewise smooth 
system, without superimposing any assumption on what kind of dynamics near E (or 
sliding motion on E) should have been taking place. The only datum for us is the original 
piecewise smooth system, and the dynamics inherited by it. 


1. Introduction 

Consider the following piecewise smooth (PWS) system: 

(1.1) x = f{x), f{x) = fi{x), X e Ri , i = 1,2,3, A , 

1991 Mathematics Subject Classification. 34A36, 65P99. 

Key words and phrases. Piecewise smooth systems, Filippov convexification, co-dimension 2 disconti¬ 
nuity manifold, Runge Kutta methods, regularization. 

This work was begun while the second author was visiting the School of Mathematics of Georgia Tech, 
whose hospitality is gratefully acknowledged, and the author would also like to acknowledge the support 
of the GNCS for the airfare. The first author also gratefully acknowledges the support provided by a Tao 
Aoqing Visiting Professorship at Jilin University, Changchun (CHINA). 

1 



2 


LUCA DIECI AND CINZIA ELIA 


for t G [0,T], and where, for i = 1,2,3,4, Ri C M** are open, disjoint and connected sets, 
and M” = System (|l.lll is subject to initial condition x(0) = xq, prescribed in one 

of the regions -Rj’s. In (ll.ip . for any i = 1,2,3,4, each /j is smooth in Ri, so that there 
is a classical solution in each region Ri, but the solution is not properly dehned on the 
boundaries of these regions. We assume that these regions are separated (locally) by an 
implicitely dehned smooth manifold S of co-dimension 2. That is, we have 

(1.2) S = {x G : /i(x) = 0, 


and for all x G S: h{x) = 


hi{x) 

/i2(x) 


Vhj(x) / 0, hj G , k > 2, j = 1,2, and 


T- '"j 

V/ii(x), V/i 2 (x), are linearly independent on (and in a neighborhood of) S. It is use¬ 
ful to think of E = Si n E 2 , where Ei = { x : hi{x) = 0}, and S 2 = { x : /i 2 (x) = 0}, are 
co-dimension 1 manifolds. Finally, without loss of generality we will henceforth use the 
following labeling of the four regions Ri, i = 1,2, 3,4: 

Ri : when hi < 0 , /12 < 0 , R 2 : when hi < 0 , /12 > 0 , 

(1.3) 


: when hi > 0 , h 2 < 0 


R 4 : when hi > 0 , h 2 > 0 


For later use, we will also adopt the notation Sj ’^2 to denote the set of points x G Ei or 
S 2 , for which we also have h 2 (x) ^ 0 or hi(x) ^ 0. E.g., = {x G Ei and h 2 (x) > 0}. 

Finally, we will denote with 

(1.4) w){x) = Vhi{x)~^fj{x) , i = l,2, j = l,2,3,4, 

the projections of the vector helds in the normal directions to the manifolds. 


For dni, a classical solution in general cannot exist on the boundaries of the given 
regions, and several concepts of generalized solution have been proposed during the years 
(see [1] for a beautiful exposition on different solutions concepts). We will restrict attention 
to Filippov solutions, m, consisting of absolutely continuous functions whose derivative 
is in the convex hull of the neighboring vector fields almost everywhere. 


1.1. Co-dimension 1. In the case of one single discontinuity manifold of co-dimension 
1, Filippov methodology has provided a widely accepted mathematical framework to un¬ 
derstand motion on the discontinuity surface. 


Consider a discontinuity manifold E = {x G : h(x) = 0 , h : M” ^ M }, separating 
two regions Ri (where h{x) < 0) and R 2 (where h{x) > 0), with respective vector fields 
/i and f 2 - Assuming that E is attractive, a condition that is satisfied when 

fi > 0 and V/ 1^/2 < 0 , x G E , 


then sliding motion on E takes place with vector field 


(1.5) 


/f := (1 - a)/i af2 , 


Vh^h 

Vh^(/i - /2) ■ 


Filippov theory provides also first order exit conditions: if one of /i or /2 (but not both) 
become tangent to E, then a = 0 or 1, E loses attractivity, and the solution trajectory 
generically will leave E tangentially (and smoothly) to enter in Ri with vector field fi, 
or in i ?2 with vector field f 2 - Furthermore, it has been understood for a long time that 
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the limiting behavior of the iterates obtained with Euler method near S leads to the 
selection of the Filippov sliding vector field itself (e.g., see [26], Chapter 3, Section 1.1) 
whenever S is attractive, and that the Euler iterates leave a neighborhood of S when S 
loses attractivity. 


1.2. Co-dimension 2. However, when S is of the form (|1.2p . an obvious lack of unique¬ 
ness (in general) arises in the construction of a Filippov vector field sliding on S. In fact, 
on S, Filippov methodology now leads to the requirement that the vector field satisfies 
the following (for positive values of Ai,..., A 4 ): 


/f := Ai/iA 2/2 + A 3 / 3 -h A 4/4 , where 



( 1 . 6 ) 


which is clearly an underdetermined system of equations. Indeed, even when S is attrac¬ 
tive, in general ()1.6p has a one parameter family of solutions and hence of possible Filippov 
sliding vector fields. 

Remarks 1.1. 

(i) Characterization of attractivity of S in the present co-dimension 2 case is consid¬ 
erably more elaborate than in the case of co-dimension 1. We will assume that S 
is either attractive by subsliding (see [TO] ) or attractive by spiralling (see 

the form recalled below in Definition \l.^ 

(ii) In the present case, the general lack of uniqueness is not resolved by considering 
the limiting behavior of the Euler iterates near S. In fact, as already noted in |15j . 
the limit of the Euler iterates selects one specific element of the one-parameter 
family of Filippov solutions. 

Definition 1.2. T, (or a portion of it) is attractive upon sliding if it is reached in finite 
time by solution trajectories for any given nearby initial condition, and further there is 
sliding motion towards S along at least one of the When there is sliding motion 

towards S along all of the £^2 then we say that S is nodally attractive. Y, (or portion 
of it) is attractive by spiralling ifY is reached in finite time by trajectories for any given 
nearby initial condition, and there is clockwise or counter clockwise motion around Y (and 
no sliding on 51^2^ fof the funetions hi{x) and /i 2 (x). 


When Y is attractive, in the literature there have been at least two systematic proposals 
to select the coefficients Aj’s in (II. 6 h . leading to sliding vector fields of Flilippov type on 
Y: the bilinear and the moments vector fields. The former has been extensively studied. 
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see [DEI HQ [IT] for example, and it consists in choosing the sliding vector field in the form 
fB ■■= (1 - «)[(! - /3)/i + /3/2] + a[(l - /3)/3 + /3/4] , 



(1 - a)(l - (3) 



rcj; W 2 wl w\ 


(1 — a)l3 


'0 

2 2 2 2 

Wi ^t?2 


a{l - fi) 


0 


af3 




The moments method has been recently introduced in |9| and consists in solving the linear 
system in (11.61) . by appending to it the extra relation 


( 1 . 8 ) 


\di —d2 —da ^ 4 ] 


Ai' 

A2 

A3 

A4 


= 0 , di = 


\Wi\\2 


i = 1,2,3,4 . 


The two choices above generally lead to distinct solution trajectories, a chief difference 
between them being the behavior of the bilinear and moments trajectories at so-called 
generic tangential exit points. These were introduced in m, and are values of x E S, 
where one (and just one) of the sliding vector fields on or is itself tangent to S 
{exit vector field). The moments vector field automatically aligns with the exit vector 
field, whereas the bilinear vector field does not. 


An important question, and one which we will try to answer in this work is: what 
should happen to trajectories of dm) in a neighborhood of S, when S loses attractivity 
(according to generic first order exit conditions)? 


1.3. Regularize. If natura non facit saltun^, then (jl.ip is either a description of an 
innatural system or a wrong model. Probably, whenever it is set forth, it is neither 
innatural nor wrong, though it may be a little of both things, perhaps because the model 
should be complemented by some missing information (the 20 years old exposition of 
Seidman in |23| is still worthwhile reading). Regardless, the above 18th century motto 
suggests considering a regularized version of (11.11) , by replacing it with a smooth differential 
system. Of course, we must assume that we do not have knowledge of where (11.11) comes 
from, if from anywhere at all, otherwise we should surely use this knowledge. However, 
in the absence of further insight, it is not obvious how one should globally regularize the 
system, and several possibilities for globally regularizing the system have been proposed 
in the literature. 

Arguably, the most studied regularization techniques are what we may call the “singular 
perturbation” and “sigmoid blending” techniques. The paper m was the first seminal 
work on the singular perturbation approach, followed by the more recent works [TB], as well 
as US], US], HZ] and, in the context of gene regulatory networks, [20]. The first systematic 
exposition of blending was the beautiful work [D. But, in the end, these techniques are all 
rather similar and amount to a regularization of the bilinear vector field. This can be done 
locally, just in a neighborhood of S, say using a cutoff function (as we do in Section [2|), 

^Linnaeus’ Philosophia Botanica, 1751 
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or more globally, perhaps through use of hyperbolic arctan functions to connect different 
vector fields. 

In this work, we are exclusively concerned with local regularizations, hence those that 
alter the given problem only in a neighborhood of S. In this case, the above mentioned 
regularization proposals share some common traits, the most important ones being that, 
outside of a neighborhood of S, the regularized vector field effectively reduces to the 
original vector fields, and that the regularization depends on a small parameter (or several 
small parameters) in such a way that as the parameter(s) go to 0, the neigborhood collapses 
onto S. The first fact is surely a reasonable property, since, away from S, there is a well 
defined smooth vector field depending on where the trajectory is, be it one of the original 
/i, 2 , 3,4 or one of the Filippov sliding vector fields on the surfaces Si_ 2 - The second fact 
may be a bit more controversial, since the regularized trajectory will typically select a 
specific sliding motion on S (when S is attractive), as the neighborhood collapses onto 
S; however, as it is well understood, and as we will also see, different regularizations do 
behave differently. As noted by Utkin ( [26] ), given that in a neighborhood of S there are 
non-unique dynamics, the inherited dynamics on S (sliding motion) does depend on the 
choice of regularization. But, this being the case, is there an appropriate way to evaluate 
different regularizations? To answer this question, we must first decide how we should 
evaluate dynamics. 

1.4. Evaluate dynamics. The first observation is that when we evaluate the dynamics 
of (jl.ip we should distinguish between the two cases (i) and (ii) below. 

(i) The PWS smooth system 0 is just a “convenient” formalism: there is a “true” 
smooth system, defined globally, but it is simpler to replace it by a PWS one. For 
example, this is the case for problems arising in gene regulatory networks (|19j. 
[5]). In these cases, one knows what is the desired behavior in a neighborhood 
of S: it is the behavior of the original problem! However, one must be careful 
in replacing the true system by (II.Ih . since it is not clear that the dynamics of 
the purely PWS system reflect the dynamics of the original problem; this was 
already noted in works such as m and [23] relative to a discontinuity surface of 
co-dimension 1. But, doesn’t this discrepancy mean that representing the original 
problem as a PWS one was not an appropriate modeling simplification in the first 
place? 

(ii) The problem arises as piecewise smooth problem, or we do not have sufficient 
knowledge of an underlying “true” problem (if any); for example, in bang-bang 
control, or in dry friction models. In our opinion, in these cases when there is no 
(knowledge of an) underlying “true” dynamics that one is trying to recover, the 
choices we make must be consistent with the PWS formulation. This is the case 
in which we are interested. 

The above consideration (ii) leads us to restrict to the model (II.ip as the system we 
are given and it is this system that we will study. This realization motivates us (and it 
has motivated us for several years) to look at the global properties of the discontinuity 
surface, and to decide on what is appropriate based on whether or not the surface attracts 
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the dynamics of the PWS system for initial conditions off the discontinuity surface itself. 
In our opinion, it is not easy to justify studying (jl.ll) by assuming a specific form of an 
underlying system from where (ll.ljl arises as some form of limiting process. In other words, 
whereas it is surely legitimate (and by no means trivial) to study the limiting behavior of a 
specific choice of regularized vector field defined in a neighborhood of S, this may actually 
have a restricted scope of applicability compared to (HH). As already noted by Utkin 
(see [25l p.44]), once one has replaced (II.Ij) with a smoothed version of it, the stability of 
sliding motion on S is inherited by that of the dynamics of the regularized field. But, is 
this consistent with the formulation (jl.ip ? 

For the given PWS system dni), we believe that it is its dynamics near the discontinuity 
manifold that determines the appropriate behavior of a trajectory. This dynamics is 
what we will try to capture in this work. Moreover, our main interest is in the case 
when the discontinuity manifold transitions from being attractive to not attractive for the 
trajectories of the original discontinuous system: in this case, will (or should) a trajectory 
(even an ideal trajectory, sliding on the manifold) feel this loss of attractivity, and hence 
leave a neighborhood of S? 

Therefore, without assuming any form of idealized motion on S, we can reformulate our 
main task as follows: 

“how should we evaluate the dynamics of dni) in a neighborhood of S ?” 

In principle, one may want to do this by studying the dynamics of a regularized problem, 
and we have already mentioned some possibilities, such as sigmoid blending and singular 
perturbation techniques. E.g., see [U [la IlSl III O [24]. We emphasize once more that 
these choices (as noted by Alexander and Seidman for blending, and by Teixeira et alia, 
for singular parturbation) remove the ambiguity of how sliding on S occurs, but these 
choices are effectively modeling assumptions, and we should ask if they render a behavior 
of the dynamics on S that is consistent with that of (II.1|) . 

Other possibilities have also been set forth in the literature, see [T] for further references. 

(a) Euler broken line approximation. This is simple to do, and it consists in replacing 
(inj by a Euler method approximation with constant stepsize, call it r. We 
have experimented extensively with this technique, see below. [The eventuality, 
of probability 0, that an Euler iterate lands exactly on a discontinuity surface can 
be handled in different ways; e.g., by randomly selecting one of the neighboring 
vector fields, or by retaining the last vector field used.] 

(b) Hysteresis (or delay) approximation. This approach appears in [26] for a surface 
of co-dimension I. For the case of S of co-dimension 2, it has been studied first 
in |2] and then in |8], always in the case of nodally attractive S. The idea here is 
that one has a region 1 /^ around S (called a chatterbox in |2]), and uses the same 
vector field, say /i, not only in the region Ri, but until the boundary of 11 ^ in a 
different region is reached; at that point, a switch to the appropriate vector field in 
the new region is performed. The rationale for this approach is that one does not 
notice immediately that a discontinuity surface is reached, but there is a “delay” 
in appreciating this fact. 
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(c) Replacing (11.11) with a stochastic DE of Ito type: x = fi{x)dt+edwt. Again we note 
that there is zero probability of landing exactly on the discontinuity surface(s). 
The interesting feature of this approach is that it is bound to sample different 
vector fields around S. The disadvantage is that it makes quantitative predictions 
possible only in a statistical sense. 

With the exception of (c), the other choices above effectively replace the original PWS 
with another deterministic dynamical system, possibly discrete (as in the case of Euler 
method). But, unfortunately, these new systems have their own dynamics, and it is unclear 
whether or not these are consistent with that of m- See Section H] for ample illustration 
of this fact. At the same time, each of the cases above has some distinguished features, 
that we will attempt to retain. 


1.5. Proposal and Plan. Our proposal to evaluate the dynamics of (11.111 in a neighbor¬ 
hood of S is to: 

“consider the Euler iterates with random steps, ” 

with steps uniformly distributed about a reference stepsize. In other words, we want to 
retain the simplicity of looking at Euler iterates, but aim at retaining a certain amount 
of randomness in the process to avoid getting trapped by the purely Euler dynamics. 
Eurthermore, in this work we will restrict to well-scaled vector fields fi, i = 1,... ,4, with 
none of the tc)’s in (II.4h exceeding 1 in absolute value. The reason for this is to avoid 
the numerical trajectory going too far away from S, and at the same time to attempt 
retaining the flavor of a hysteretical trajectory. 

We will complement our experiments made with the above strategy, by also using other 
approaches. For example, Euler method with constant stepsizes, singular perturbation 
regularizations, and also numerical integration of o performed with variable stepsize 
integrators. 

We emphasize that our goal is to reach some conclusions insofar as what should happen 
to a solution trajectory of (jl.ll) in a (small) neighborhood of S. We are not compar¬ 
ing methods, or different recipes of ideal sliding motion, but simply trying to evaluate 
the dynamics of (HI) in the most plausible and honest way we can think of, without 
superimposing on o any extra modeling assumption. 

Finally, one more caveat. Our examples are all of systems in and M^, and the 
discontinuity surfaces Si and S 2 are planes; this makes it easier to visualize and understand 
things. Higher dimensional state space, and non-planar discontinuity surfaces, can surely 
bring new phenomena into play, but we have no reason to suspect that the basic picture 
that emerges in our study, with the dichotomy between attractivity and lack of attractivity 
of S, will be modified substantially. 

The remainder of this work is structured as follows. In Section [21 we consider a pro¬ 
totypical regularization of the bilinear vector field (jl.7p . and give sufficient (and sharp) 
conditions guaranteeing that the regularized solution converges to a sliding solution on S 
according to (II.7p . In Section [3l we give some details of how we implemented the above 
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mentioned proposal, particularly of what practical criteria we adopted to detect “exit¬ 
ing” from a neighborhood of S. In Section [U we illustrate through several examples the 
different things that can happen, and their dependence on the adopted simulation choice. 


2. Bilinear interpolant regularization: solutions behavior and fast slow 

DYNAMICS 

Space regularizations are often employed in literature as an analytical mean to model the 
switching mechanism of a discontinuous system, see [21 |T3l UHl |T71 HSl [HI [25] . Typically, 
these regularizations are one parameter families of vector fields with different time scales 
in a neighborhood of the sliding surface S, namely a slow dynamics tangent to S and a fast 
one normal to S. In what follows, we consider regularizations for Filippov discontinuous 
systems as in (inD, but other approaches are available in the literature for non-smooth 
systems that are not of Filippov’s type, for example control systems with nonlinear control, 
as in [25l [26|. When the regularization parameter goes to zero, the regularized solution 
converges to a solution of Filippov’s differential inclusion ( [15[ Theorem 1, §8]). It follows 
that, if S is a codimension 1 discontinuity surface, for any regularization that satisfies 
the assumptions of m Theorem 1, §8], the corresponding solution will converge to the 
unique sliding Filippov solution on S as the regularization parameter goes to zero. But 
if S has codimension k > 2, then the ambiguity of Filippov’s selection will reflect also in 
the amibiguity of the limit of regularized solutions (in other words, the limit will depend 
on the chosen regularization). 

Our goal in this section is to study the limiting behavior of the solutions of a certain 
regularization, namely the bilinear interpolant, or simply bilinear, regularization (|2.2I) 
below. This regularization (or a close relative) has often been employed and studied in 
the literature; see [Il[2l[22l|23lll3llI71[l6l[l8]. More specifically, we will study when the 
solution of the bilinear regularization converges to the solution of a particular selection 
of the Filippov’s vector field: the sliding bilinear vector field (|1.7p . When S is nodally 
attractive, this convergence has been shown in [T] Theorem 5.1], and -under the same 
assumption- in m it is shown that the bilinear regularized vector field converges to the 
bilinear sliding vector field. In what follows, we relax the hypothesis on attractivity of 

S, and simply assume that S is attractive in finite time, i.e. all trajectories with initial 
conditions in a neighborhood of S will reach S in finite time. This can be achieved either 
upon sliding along Si and/or S 2 (S is attractive upon sliding), or spiralling around S (S is 
spirally attractive); see Definition [T21 For later reference, and under the stated attractivity 
assumptions of S, we note that the algebraic system (HZj) has a unique solution {a, f3) in 
(0,1)2. 

For simplicit30, we assume that S is the intersection of the hyperplanes Si = {x G 
M”’]®! = 0} and S 2 = {x G M"'|x 2 = 0}. Then, we consider the e-neighbor hood S of 

T, : S = {xi,X 2 : —e < xi,X 2 < e}, and two smooth (at least C^) monotone functions 
(the functions of course depend on e, but for notational simplicity this dependence on e is 


fact that can be always locally guaranteed 
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omitted) a, (3 : M —>■ [—1,1] interpolating at ±e, as follows: 

(2.1) a{e) = /3(e) = 1, a(—e) = f3{—e) = 0, a'{xi), f3'{X 2 ) >0, for — e < xi,X 2 < e . 


We call bilinear regularization the following one parameter family of vector fields 


feix) = (1 - a(a;i)) [(1 - /3 (x2))/i(x) + P{x2)f2{x)] + 
a(xi)([(l - I3{x2))h{x) + /3(x2)/4(x)] . 


To be specific, in what follows, we have taken 
(2.3) 


q;(xi) 


I+ 1^(3-(f) 


xi > e 
xi G [-e,e] 
xi < —e 


I3ix2) 


+ f(3-(f)^) 


X 2 > e 
X2 G [-e,e], 
X 2 < -e 


but other choices of a monotone interpolant could be considered. Clearly, a choice of 
two different parameters for a and /3, respectively Cq, and e^, could also be considered (e.g., 
see [22]), and this would be justified if, for example, it is known a priori that the trajectories 
of un underlying physical system approach the two surfaces Si and S 2 at different rates. 
Naturally, the choice of different parameters might lead to different qualitative behavior 
of the corresponding solutions, as we will see in Section jH 


Remark 2.1. Often (e.g., see mEaEsi;, a simpler regularization is adopted, namely 
q;(xi) = , and similarly for /3(x2), possibly with different and ep. Provosition \2. (j\ is 

not meaningful in this case, whereas Proposition \2. 7| holds with essentially the same proof. 
See also Remark \2.4\ 


Our goal is to study the behavior of solutions of x = /g(x) for x £ S, and to see when, 
for e —)• 0, they converge to the sliding solution on S with vector field fs as given in (ll.7p . 
In order to do so, we split x = /g(x) into fast and slow motions. For (xi, X 2 ) in S, we can 
rewrite (|2.2I) with respect to the variables (a, /3, X3,. .., x^). For the sake of more compact 
notation, we let y = (xs,... ,Xn) and further y = g{a,l3,y). In these variables, the full 
system rewrites as 

a = ^e{f^{a,(3,y) := ^gi{a,l3,y) 
dx\ dxi 

(2-4) o d/3 d/3 

y = 9{a,/3,y) 

where e, is the standard i-th unit vector in i = 1,2. Notice that and depend 

2 

on e and are strictly positive (inside S); from (12.31) . they are equal to ^(1 — ^), with 
i = 1, 2. We refer to a and /3 as the fast variables and to y as the slow variable. We denote 
the solution of (|2.4I1 as (ae(-))/3e(‘)) ye('))- 

Now, using (|2.3p . in (|2.4p . because of monotonicity of a and (3, we can rewrite xi as a 
function of a and X 2 as a function of /3. From ()2.3p . let z = ^ and rewrite the first of (|2.3p . 
for xi G [— e, e], as: — 3z + 4 q: — 2 = 0. Using Vieta’s substitution z = tc + ^, we get the 
equation re® + (4a —2)t(;^ + l = 0. Then = (1 — 2a) ± 2Wa — o? and its third roots have 
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modulus one, so that z = w-\-^ = w + ^ is real. Let (p{a) = Arg ^(1 — 2a) + 2 iy/a — a^j , 

• / ipi^oc) 4 

then for a in [ 0 , 1 ], </?(«) is in [ 0 ,vr], and for w{a) = e*''“3 ^ 3 ^), the corresponding z{a) 
can be rewritten as z{a) = 2 cos(^^^ + | 7 r) and satisfies z( 0 ) = —1 and z{l) = 1 , this 
z(a) is the function we are looking for. Same reasoning applies for j3. Then system (|2.4I) 
rewrites as 


(2.5) 


ed = ^ (l - 4 cos^ /3, y) 

e/5 = ^ (l -dcos^ ^2Ma,l3,y), 

y = 9{(^,P,y) 


with ip{'y) = Arg ^(1 — 27 ) + , Notice that the function ^1 — 4cos^ + 1”^)) 

is strictly positive for 7 G (0,1) and it is 0 at 7 = 0,1. Following standard approaches for 
singularly perturbed systems, we set e = 0 in ()2.5I) and obtain the following system 

0 = gi{a,j3,y) 

(2.6) 0 = 52(a,/3,y) 

y = 9{a,P,y) 


Notice that solutions of (12.61) are sliding solutions on S with bilinear vector field fs as in 
dEZI). Let (a*(y), (3*(y)) denote the solution of 


gi{a*{y),/3*{y),y) =0 

g 2 {a*{y),/ 3 *{y),y)) =0’ 


(Recall that, under the assumptions of attractivity by sliding or by spiralling, {a*(y), /3*{y)) 
is the unique solution of (j2.7|) in [0,1]^). 

Our goal in this section is twofold: to see if solutions of (|2.4[) converge to solutions of 
(j 2 . 6 [) as e —)• 0 and to explore the behavior of solutions of ( 12 .4p in the neighborhood of 
generic exit points. 


2.1. Asymptotic behavior of the regularized problem for S attractive in finite 
time. From ()2.4n . we introduce the time variable t = || and consider the fast system 


a' = - 4 cos^ ) 91 (a, / 5 , y) 

/5' = (^1 - 4cos2 92 ( 01 , P,y), 


where the “prime” denotes differentiation with respect to r, and y is considered as a vector 
of parameters. Notice that the solution (a*(y), f3*(y)) of (12.7p is an equilibrium of (12.811 . 
We denote with (a(t,y), I3{t,y)) the solution of (|2.8p with initial condition (ooj/^o)- 
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Remark 2.2. In case we take different and ep in (j2.3p . we can write ep = rjCa, perform 
the change of time variable r = and obtain (similarly to ()2.8p ) the fast system 


(2.9) 


/ 3 ' = 


a' = ( 1 — 4 cos^ 
1 


+ ^^ ) ) giia,P,y) 


V 


1 - 4cos2 ) 92ia,/3,y) • 


The following result by Artstein is a stronger version of a classical result in singular 
perturbation theory (O Theorem 2.1]). 

Theorem 2.3. Assume that 

(i) the solution {a*(y),/3*(y)) of (12.7p is continuous in y; 

(ii) (a*(y),/3*(y)) is a locally asymptotically stable equilibrium of the fast system (|2.8I) .' 

(hi) the initial condition (oQ) /^O) yo) of (12.4p is such that the co-limit set of{a{t, yo), /3{t, yo)), 

is (a*(yo),/3*(yo)); 

(iv) the problem y = g{a*{y), I3*{y),y), y{0) = yo, has a unique solution, denote it as 

yo{t)- 

Then, the solution of (12.4p with initial condition {oq, l3o,yo), is such that as e ^ 0.- 

a) y{-) converges to yo{-), uniformly in time on intervals of the form [0,T]; 

b) (a(t),/3(t)) converges to {a* {yo{t)), {yo{f))) uniformly in time on intervals of the 

form [(5,T], (5 > 0. 

□ 


Therefore, when conditions (i)-(iv) of Theorem 12.31 are verified, the solution of the 
regular system converges to the sliding solution with bilinear vector field ()1.7p . 


Remark 2.4. If we used (see Remark \2. 1\) the regularization a{xi) 
and let T = ■^, then we would have obtained the system 


( 2 . 10 ) 


a’ = gi{a,(3,y) 
13' = g2{a,f,y) 


Xi+t 
2e ’ 


P{X2) = 


instead of ()2.8p. Now. (]2.10p is precisely the “dummy” system, of m and [H] . How¬ 
ever, (12.8|) is not orbitally equivalent to (|2.10p . As we will see in Proposition |g. 7[ under 
appropriate conditions of attractivity ofT,, {a*(y), (I* (y)) is an asymptotically stable equi¬ 
librium for both (j2.8p and (I2.10h . However, condition (iU) in Theorem, \2.A implies that 
the limiting behavior of the solution of (|2.4p . depends also on the basin of attraction of 
{a*(y), ft*(y)) and this, in general, is not the same for the two systems. As an illustration 
of this, see Example \4-f^ in Section System (|2.8p is the system we need to study in 
order to understand the limiting behavior of the regularized solution in the case of the 
regularization dZH). Moreover, we note that if we had used the regularization with 
different parameters Cq, and eg, namely a{xi) = and (3{x2) = then we would 

obtain a system not orbitally equivalent to (|2.inp . 
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In what follows, we first assume that S is attractive in finite time (upon sliding or 
spirally) and we want to verify if and when (i), (ii) and (hi) are satisfied. This in turn 
will imply that solutions of the regularized problem converge to the sliding solution on S 
with vector field Jb as given in (ll.7jl . The following hold, both when S is attractive upon 
sliding or spirally. 

(a) System (|2.8jl has a unique equilibrium (a*(y),/!*(?/)) in (0,1)^. 

(b) The functions a* (y), (3* (y) are smooth functions of y. (For S attractive upon 
sliding, this is in m Theorem 8], and the case of S spirally attractive is analogous.) 
Notice that this point (b) implies (iv) of Theorem 12.81 

Remark 2.5. We emphasize that the assumption that S is attractive in finite time is 
sufficient but not necessary for the uniqueness of {a*{y), fi*(y)); e.g., see [TO]. As we will 
see in Section^^ this in turn will impact the behavior of the regularized solution, that might 
remain close to S even when S is not attractive. 

In Proposition 12.61 and Proposition l2.71 we study the dynamics of ()2.8p to verify asymp¬ 
totic stability of 

Proposition 2.6. The phase space of (12.8p is the square Q = [0,1]^. Moreover, the 
boundary of Q, denoted as dQ, is an invariant set for all values ofy. The vertices of Q 
are equilibria. If any of the sliding vector fields ^ defined (i.e., if there is sliding 

on any of'Ef 2 )> Ihe corresponding value of {a{y), I3{y)) in (|I.7p is an equilibrium of (12.8p 
for y = y. 

Proof. For a = 0,1, or /3 = 0,1, it must be xi = — e, e or X 2 = —e, e, and hence a' = 0 
or fi' = 0. Then, for all y, the boundary of Q is invariant under the flow of (12.81) . It also 
follows that the vertices of Q are equilibria of (12.8p . Notice that the corresponding values 
of gi and g 2 at the equilibria are (rc|,rcf) for (0,0), for (1,0), (tC 2 ,rc|) for (0,1), 

and {wl,w1) for (1,1), where the w^s are defined in ()1.7I) . 

In addition to the vertices of Q, other equilibria of (j2.8l) might belong to dQ, as follows. 

• If, for y = y, there is sliding on then there exists I3~{y) such that 52 ( 0 , /3~ iy),y) = 
0 and (0,/3“(y)) is an equilibrium of (12.8h . The vector field (II.7h for {a,/3) = 
(0,/3“(y)) is 4-(y). 

• If, for y = y, there is sliding on then there exists I3~^{y) such that 52 ( 1 , /3~^{y),y) = 
0 and (1,/I'''(y)) is an equilibrium of (|2.8[) . The vector field (jl.71) for (a,/3) = 
(0,/3+(y)) is 4+(y). 

• If, for y = y, there is sliding on for y = y, then there exists a~{y) such that 
gi{a~{y),0,y) = 0 and (a“(y),0) is an equilibrium of (|2.8jl . The vector field (II.7h 
for {a, 13) = (a"(y),0) is f^-{y). 

• If, for y = y, there is sliding on for y = y, then there exists a^[y) such that 

g 2 {ot'^{y), l,y) = 0 and 1) is an equilibrium of (|2.8I) . The vector field ()1.7[) 

for (a,/3) = (a+(y),l) is /s+(y)- 

□ 
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Table 1. 



i = 1 

i = 2 

i = 3 

i = A 

wl 
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Figure 1. Proposition 12.71 Reference configuration under which S is at¬ 
tractive in finite time. 
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! y 

X 
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In Proposition 12 . 71 we give sufficient conditions for (a*,/3*) to be locally asymptotically 
stable. 

Proposition 2.7. Assume that S is attractive in finite time upon sliding along at least 
two of the ihat there is attractive sliding along these codimension 1 surfaces. 

Then {a* (y), f3*(y)) is exponentially asymptotically stable for (12.81) . 


Proof. We prove the result for a particular configuration of vector fields in a neighborhood 
of S. The proof for all the other configurations is analogous. We consider the case in Figure 
[H where S is attractive upon sliding along ^ characterized by the signs of the wys 
in Table [H where the wys are dehned in (II.4h . 

S attractive implies that (12.41) has a unique equilibrium (a*,/3*) E (0,1)^. To study the 
stability of (q;*,/ 3*) we consider the Jacobian matrix of (j2.8p . evaluated at (a*,/3*), that 
is we look at the eigenvalues of 


( 2 . 11 ) 

/(l-4cos2(*!>+|,))^(a-,r) (l- 

^l_4co82(£m + |^))^(a-,/C) (l- 

/l — 4cos^(^^)p + jTt) 0 

I 0 1 — 4cos^(ll(^ + Itt) 


4co52(^ + l>r))^(a-,r)j 

I 7 7-= 
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Then, as in the proof of m Theorem 8], it follows that det(J) > 0, and hence det(J) > 0. 
We now show that trace(J) is negative. We write explicitly the diagonal elements of J 

= -(1 - - /3*wi + (1 - 

/?*) = —(1 — a*)wi + (1 — a*)w 2 — a*w‘^ + a*w\. 

and the equilibrium (a*, = ( ci(/ 3 ^)-C 2 (r) ’ C'i(/3*) = (l-^*)ini + 

I3*W2, C2{/3*) = {l—/3*)wl+/3*wl, Li{a*) = {l—a*)wi+a*w‘^, L 2 {a*) = {l—a*)w 2 +a*w‘l. 
Notice that 0 < a* < 1 together with C'i(/3*) > 0 (see Tabled]) implies that ^(a*, /3*) < 0. 
Similarly, 0 < /3* < 1 together with ^ 2 ( 0 *) < 0 implies ^(a*,/3*) < 0. This together 

with — 4cos^(^^^ + Ivr)^ > 0, gives the sought result. □ 

The assumptions in Propositions 12.71 exclude the case of S spirally attractive and S 
attractive upon sliding on just one of the For each of these cases, examples can 

be given where the equilibrium (a*,/3*) of (12.8p is unstable even if S is attractive. See 
Example 14.31 in Section d] where, with e^g = lOco, the equilibrium {a*, j3*) of the fast 
system is unstable while S is attractive. 

2.2. Behavior of the regularized solution in the neighborhood of exit points. 

Here, we are interested in studying how solutions of (12.4p behave when S looses attractiv- 
ity. Assume first that S is attractive upon sliding. We will focus on the case when S looses 
attractivity at so-called potential tangential exit points and at potential non-tangential exit 
points, defined next. 


Definition 2.8. Let x be a point on S. We say that x is a potential tangential exit point 
if, at X, one of the vector fields ^ 2 to S (and hence it belongs to the Filippov 

convex combination). We say that x is a potential non-tangential exit point if, at x, one 
of the vector fields fj{x) is tangent either to Si or to S 2 and it points away from S. 

Simplification: S is a curve. For simplicity, below we restrict to the case n = 3, so 
that, in our case of S = {xi = X2 = 0}, S is a (portion of the) xs-axis. (Of course, in 
general, when the co-dimension 2 discontinuity manifold S is immersed in M”, with n > 3, 
we expect that exit points themselves will lie on a (union of disjoint) (n — 3)-dimensional 
manifold(s) of codimension 3. With the understanding that exit points are now lying on 
these manifolds, our description below still qualitatively holds). Also, we consider only 
first order phenomena i.e., only potential tangential (respectively, non-tangential) exit 
points X, such that the derivative with respect to time (evaluated along the trajectory) of 
/si 2 ^^) (respectively, fj{x), j = 1,..., 4) is different from zero. As a consequence, for us 
exit points are isolated points on S. In this case of n = 3, and insofar as the behavior of 
the regularized solution in the neighborhood of tangential and non-tangential exit points, 
the following phenomena can arise. 
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2.2.1. Tangential exit points. What occurs in this case depends on the number of equilibria 
of ( 1 ^. 

Suppose that the regularized solution enters the neighborhood of a tangential exit point 
X = ( 0 , 0 , X 3 ), and without loss of generality let fj.+ {x) be the vector field tangent to S, 
and let S be attractive for X 3 < X 3 . Then either (a) or ( 6 ) below can happen. 

(a) {a* {X3), (3* (x^)) = {a~^{x 3 ), 1) and the regularized solution exits the neighborhood 
of S to enter a neighborhood of ; 

(b) (a*(a: 3 ),/ 3 *(x 3 )) / (a"''(x 3 ),l) and these are the only solutions of (12.71) for y = 

X 3 = X 3 ', a first order analysis guarantees that the solution that at X 3 is equal 
to 1) is in the interior of [0,1]^ for X3 in a right neighborhood of X3. 

Nonetheless, {a*{x 3 {t)), I 3 *{x 3 {t))) retains its asymptotic stability. It follows that 
the regularized solution remains close to S and for e —)• 0 it will converge to the 
solution of the bilinear vector field on S, still well defined even though S is not 
locally attractive anymore. Eventually the regularized solution x{t) will leave a 
neighborhood of S if one of the two following phenomena takes place: 

(i) the equilibrium of (12.81) reaches a fold bifurcation value; 

(ii) {a*{x{t)), I3*{x{t))) is on the boundary of Q foi t = t, and this can happen 
only if one of the vector fields ^(x(t)) is tangent to S, say 

In case (ii), the regularized solution will leave the neighborhood of S to enter a 
neighborhood of But, in case (i) it is not generally possible to predict how 
the regularized solution will leave a neighborhood of S. 

The phenomena (a) and (b) above do not depend on ^ and ^ in (j2.8l) . but only on the 
solution of the algebraic part in ()2.6p . Hence any choice of functions a and /3 in (12.2j) that 
satisfies ( 12 .Ij) will lead to the same phenomena. 


2 . 2 . 2 . Non tangential exit points. Suppose that the regularized solution enters the neigh¬ 
borhood of a non-tangential exit point, x = (0,0, X 3 ). In this case, although there is only 
one solution of (| 2 . 6 p . different things can happen depending on the functions and 
in ( 12 . 8 p . 

Suppose that /i is the vector field verifying the conditions for non tangential exit. 
Then (q;*(x 3 ), / 3 *(x 3 )) is still the only equilibrium of (|2.8I) in the interior of [0,1]^ and it is 
asymptotically stable. However, the equilibrium (0,0) undergoes a bifurcation at X 3 = X 3 
and, for X 3 > X 3 , there is a neighborhood of (0,0) in Q such that all solutions in that 
neighborhood are attracted to ( 0 , 0 ). 

As far as the regularized solution, (i) or (ii) below may occur. 

(i) The regularized solution x{t) might remain close to S, in agreement with Theorem 

EH 

(ii) The regularized solution will leave a neighborhood of S to enter i?i, if the corre¬ 
sponding solution of ( 12 . 8 h enters the neighborhood of solutions that reach ( 0 , 0 ). 
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In other words, the behavior of the regularized solution depends on the possible bifurca¬ 
tions of (a*,/3*) and of (0,0). Moreover, we can choose the functions a and /? so that the 
terms and will make either one of (i) or (ii) above take place. See Example 14.21 in 
Section 3| for an illustration of this fact. 

2.2.3. Exit point in the case of spiral dynamics. We consider the case of S spirally attrac¬ 
tive and, based on the results in [7], we propose the following definition of potential spiral 
exit point. 

Definition 2.9. Let x & T, and assume that there is spiral dynamics around S. We say 
that X is a potential spiral exit point if x is such that 

w‘f{x)wl{x)w‘l{x)wl{x) ^ 
rcj; (x)w^ (x)wj (x)w 2 (x) 

Again, we are only concerned with first order phenomena, i.e., phenomena such that 
the derivative with respect to time of the quantity on the left hand side in Definition 12.91 
evaluated at x(t) = x, is different from zero. 

The equilibrium (a*,/3*) of (|2.8I) is always unique and well defined in (0,1)^ for all 
X 3 E M as long as there is spiral dynamics around S. The attractivity of S though, 
does not imply the stability of the equilibrium. On the other hand, (a*,/3*) might be 
stable when S is not attractive. As an illustration of both phenomena, see Example 14.31 
Furthermore, when the equilibrium of ()2.8|) is unstable and the trajectories of (j2.8p reach 
(possibly in finite time) one of the equilibria on the boundary of Q, it is the local dynamics 
around S that determines the behavior of the regularized solution. This is well illustrated 
in Example 14.31 Figure [TOl for = lOea. 

Based upon the above described situation insofar as different behaviors of the regularized 
solution, it is natural to ask how should the solution of the original discontinuous system 
O) behave once a potential first order exit point is reached, and how we should infer 
this. We discuss this aspect in the next section. 

3. Numerical simulations for the discontinuous system 

Frequently, Euler’s method has been used as a mean to approximate the behavior of 
solutions of (jl.ip in a neighborhood of the discontinuity surface. In [15] (proof of Theorem 
1, page 77), Filippov showed that the solutions obtained with Euler’s method converge to 
one of the solutions of Filippov’s inclusion when the discretization stepsize goes to zero. In 
[I] , Alexander and Seidman -in the case of a nodally attractive codimension 2 discontinuity 
surface E- consider a chattering trajectory that evolves in an e-neighbor hood of S. The 
trajectory x^ is obtained by considering the Euler approximation of the solution, and the 
fact that S is nodally attractive guarantees that -for a sufficiently small stepsize r- the 
Euler’s approximations remain in an e-neighborhood of S. Alexander and Seidman in 
[T] show that every possible Filippov sliding vector field on S is realizable; i.e., given a 
solution x{t) of Filippov’s differential inclusion, there exists a chattering trajectory x^ that 
converges uniformly in t to x in a given time interval. 
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3.1. Discretization. We also use Euler’s method as a tool to understand how the solu¬ 
tions of (jl.lj) should behave in a neighborhood of the codimension 2 discontinuity surface 
S. However, we propose to do this in a new way, which enables us to make (statistical) 
inferences on the expected behavior of a trajectory. 


First, we evolve by Euler’s method several initial conditions in a neighborhood of S. 
Second, in order to mimic the non ideal behavior of a physical system, at each step k, we 
use a random stepsize such that t < < 2r, where r is a fixed reference value. In case 

one of the computed approximations falls on a discontinuity surfac^, we take a random 
perturbation of size of machine precision, so that the perturbed value belongs to one of 
the regions i?j’s. Therefore, for each initial condition xq, chosen in a r-neighborhood of 
S, we generate the following approximate solution 


(3.1) 


^fc + 1 — 


f Xk+TkfjiXk), 

\xk + dk + Tkfj{xk + 6k), 


Xk ^ , 

Xk € E, Ei^2) 


where 5^ > 0 is a random perturbation the same size of machine precision. We perform 
two types of experiments with the scheme (|3.1jl . 


(a) We generate 100 (or more) random initial conditions in a r-neighborhood of E, 
and evolve each of them according to ()3.ip on a given time interval of interest. We 
further monitor, see below, exit points for each trajectory, and perform statistics 
on these. 

(b) For a given (randomly generated) initial condition in a r-neighborhood of E, we 
generate 100 trajectories according to (13.11) . on a given time interval of interest. 
Of course, becase of the randomnes in the stepsize selection process, these will give 
approximations at different times. Thus, we further interpolate linearly the given 
trajectories on a fixed temporal grid with spacing r, that is at times 0, r, 2 r,... . 
Finally, we compute an average trajectory by averaging the obtained 100 approxi¬ 
mations on the fixed grid. We also compute exit points, etc., with respect to this 
average trajectory. 


3.2. Construction of the experiments. As previously remarked, our experiments are 
all in or and with discontinuity surfaces given by hi{x) = xi and h 2 {x) = X 2 - 
Moreover, we choose the vector fields so that: 

• All but one of the w^s are constant and they are in absolute value less than 1. 
The non-constant w) is a function of the slow variables x^,... ,Xn and it is chosen 
in such a way that E changes from locally attractive in finite time to non attractive; 

• No vector field in the Filippov differential inclusion has equilibria (or more com¬ 
plicated invariant sets) on S. 

3.3. Exits. An important aspect of our simulations will be to monitor if a Euler ap¬ 
proximation {xfcjfcgis} of (|l.ll) . computed as in (|3.1I) . leaves a neighborhood of E when S 

case that has never occurred in the several thousands experiments we have performed 
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loses attractivity. To perform this task, we reasoned as follows. (The choices below are 
appropriate for the vector fields we chose, see Section ESI). 

(a) Exit on a codimension 1 surface [Tangential Exits]. Suppose that S loses attrac¬ 

tivity at a tangential exit point for which fj.+ becomes tangent to S. We monitor 
the function /i 2 and declare that the numerical solution leaves a neighborhood of S 
by checking that h 2 {xk) > lOr. When this is the case, we define as “exit point” the 
value X = where the index k is the first index k for which /i 2 (x^) > 1.5r 

together with h 2 {xj) > r for all j > k. Note that we do this both for the 100 
trajectories generated by (El]) with different initial condition, as well as for the 
average trajectory. 

(b) Exit in one of the Rj’s [Non Tangential Exits]. Suppose, for example, that the 
non-tangential exit is into the region Ri. We observe that, at the non-tangential 
exit point, /i is tangent either to Si or to S 2 , while pointing away from S. This 
implies the existence of repulsive siding along Si or S 2 with sliding vector field /i 
at the non-tangential exit point. Hence, we still need to detect first the exit on a 
codimension 1 sliding surface, which we do as in (a). 

(c) Exit from spiral case. As long as S is spirally attractive, motion around S now 
repeatedly takes the trajectory inside all of the regions Ri,..., R 4 ^. To decide if the 
trajectory leaves a neighborhood of S (when the spiral motion ceases to be attrac¬ 
tive), we monitor the 2-norm of the vector (/ii, / 12 ) for the last computed numerical 
value in R 4 , before the solution enters the neighboring region R 2 (clockwise mo¬ 
tions around S) or R^ (counterclockwise motion). We do this for the trajectories 
associated to different initial conditions, and declare an exit when the 2-norm of 
(/ii, / 12 ) becomes strictly monotone increasing. 

3.4. Summary and limitations. As we report in the next section, based upon our ex¬ 
periments with Euler method with random stepsizes, we can thus summarize our findings. 

• All computed solutions remain in a 0(r)-neighborhood of S as long as S is attrac¬ 
tive. This confirms that, on an ideal system, sliding should be taking place along 

S. 

• All computed solutions move away from S when they reach a sufficiently small 
neighborhood of a potential exit point. 

There are important caveats to our results. 

• If the discontinuous dynamical system is the “idealization” of a known smooth 
system, the approach described above may be misleading if we want to understand 
what happens to the solutions of the original system. Eor example, if the original 
system is given by ()2.2h . as the regularization parameter goes to zero, one may 
obtain a limiting behavior which differs from the behavior of the discontinuous 
system. What this means is that in this case the discontinuous model (II.Ih is not 
a good model; as an example of this, see m- 

• Another class of systems that does not fit our analysis is given by control sys¬ 
tems with non linear controls for which the Eilippov convex combination would 
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not contain all the neighboring values of the vector fields, and hence our Euler’s 
approximations of (jl.ip are not sufficient to understand the behavior of solutions 
(see [26] for a general theory). 

4. Numerical experiments 

This section presents results of numerical simulations on several examples in and 

On each example, we compare the results of different experiments carried out with 
some of the following strategies. 

1) Random Euler. This refers to the approximation of (II.ip computed with Euler 
method with random stepsizes as discussed in Section [3l 

2) Deterministic Euler. This is the fixed stepsize implementation of Euler method 
directly on (|l.ll) (see item (a) at p. [6|). 

3) Regularized integration. This refers to the approximation of the bilinear regularized 
vector held as proposed in Section [J] by integrating the regularized system with the 
Matlab built in functions ode45, and/orode23s, and/or odel5s0 All of these are 
variable stepsize integrators, where local error control is enforced by a combination 
of relative and absolute error tolerances (Reltol and Abstol). (It is surely possible 
that different solvers may perform somewhat differently than those we have used, 
but we have no reason to suspect that a different solver will alter the outcome of 
our observations below.) 

4) Unregularized integration. This refers to the approximation of (jl.ip computed (as 
a discontinuous system) with the Matlab built in functions ode45 and/or ode23s 
and/or odel5s. 

The most noteworthy aspects conhrmed by our numerical experiments are the follow¬ 
ing. The numerical solution computed with “Random Euler” remains close to S as long 
as S is locally attractive, and instead leaves any small neighborhood of S when S looses 
attractivity. On the other hand, the approximations computed by the “Regularized inte¬ 
gration” might remain in a neighborhood of S even when S is not attractive in agreement 
with the dynamics of ()2.8p . Also, “Deterministic Euler” is not a foolproof option, since it 
superimposes its own dynamics to the true dynamics of the underlying problem. Einally, 
“Unregularized integration” may just fail when requiring stringent values of the error tol¬ 
erances, and also occasionally producing totally erroneous approximations, and cannot be 
taken as a trustworthy indicator of what should happen in our context. 

Eor our experiments, we will adopt the simplifications discussed in Section [3l As already 
anticipated in Section [3l the examples are in and the discontinuity surfaces are 
Sj = {x G M”, Xj = 0}, j = 1,2 and the rc^’s are all constants except for one of them 
that depends on x G S. Moreover, with the exception of Example 14.41 in all the examples 
below the Eilippov vector field is uniquely defined. This is not a restriction for our scopes, 

^ode45 is a Runge-Kutta integrator, suitable for non-stiff problems, ode23s is a second order integrator 
based on Rosenbrok formulas, better suited for stiff problems, and odelSs is a variable order (1 through 
5) integrator based on the backward-differentiation-formulas (BDF), also suited for stiff problems. 
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since the behavior of solutions (whether or not for the regularized system) at potential 
exit points remains ambiguous, as we will see below. 


Example 4.1 (S looses attractivity through a potential tangential exit point). Consider 
the following vector fields 


(4.1) 



( i ^ 


( 1 \ 

fl = 

0 0.9 ™2 

^ 4 X3 X4 

1 /2 — 

- 0.3 

2 , 


-^X3 + Xi 


-50:3 + 0:4 




V-fx 3 + IxiJ 


( ^ 


/ -0.25 \ 

0.9 

II 

-0.15 

-|X 3 + X4 

-|X 3 + Xi 

V-|X3 + Ixi) 


\-lx 3 + IxiJ 


S is the (x 3 ,X 4 ) plane and on it there is the unique vector field x 


( 0 \ 
0 

-|X3 + X4 

V-fx3 + |X4/ 


The 


circle 7 = {(x 3 ,X 4 ) : = 2 } is a curve of potential tangential exit points on E. 

The region inside 7 is attractive upon sliding. It is attractive upon sliding along ^,^2 
x'^+xl < 2 — ^ and it is attractive upon sliding along Sf and for 2 — ^ < < 2 , 

{see [m for theoretical studies of an attractive co-dimension 2 surface S when one of the 
’s is zero). Outside 7 , fj.- points away from T, so that E is not attractive. All solutions 
with initial condition inside 7 will eventually meet 7 . We expect trajectories of (HH) 
to leave E when they reach 7 and to start sliding along E^ in direction opposite to E. 
This behavior is well illustrated in the right plot of Figure [3 where we depict the average 
solution computed with “Random Euler” with r = 10“®. Next, we report on four more 
experiments. 


a) Random Euler. We computed the average exit point for an ensamble of 100 initial 
conditions with xi and X 2 uniformly distributed in [—t, r]^ and X 3 + X 4 = 1.7. We 
compute the exit point for every solution as described in Section\^ For r = 10“®, 
we obtain the average exit point x ~ 2.0865 with standard deviation ~ 0.0202. 
In the left plot in Figure [H we plot the average trajectory obtained with stepsize 
T = 10“®. The exit point for the average trajectory is x 2.08. 

b) Regularized Integration. For the regularized vector field (12.2p . we take a and (3 as 

in (ESI). The corresponding fast system (12.81) has a unique asymptotically stable 
equilibrium in (0,1)^, a stable node, up to x\ + x\ = p c:z. 2.225. The curve 
x\ + x\ = p is a curve of saddle node bifurcation values for the fast system and 
for X 3 + > p there are no equilibria in (0,1)^. (For this example, the behavior 

of (12. 8 h does not change by choosing different functions a and /3 in (12.3h . as long 
as they satisfy conditions m)- The behavior of the solution of the regularized 
system is well illustrated in the right plot in FigurelB The solution of the regular 
problem is computed with the Matlab function ode23s and RelTol=AbsTol=10“®. 
The continuous line in the plot is the first component of the solution, while the 
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Figure 2. Example 14.11 Left: Average trajectory obtained with Random 
Euler and r = 10“®. Right: Regularized Integration, with ei = 62 = 10“® 
and initial condition [ 10 “^ 10 “^ 0 \/L5]. 




dashed line is the second component. They are both plotted in function of z = 

c) We also computed the solution of dni) with “Deterministic Euler” and fixed step- 

size T = 10“^, and initial condition (10“^, 0, \/L7). First, we get rid of the 

transient and then plot the first two eomponents of the approximation in Figure 0 
None of the elements of the sequence generated by the forward Euler’s map is on 
El or Ti 2 and hence the computed solution is always in one of the Rj’s and the 
seleetion of the vector field is straightforward. As it is clear from the plot, we do 
not recover the dynamies of the original system. The periodic orbit in the {xi,X 2 ) 
plane is spurious and is generated by forward Euler’s own dynamics. The periodic 
orbit survives past the eurve of potential tangential exit points and indeed the third 
and fourth eomponent of the plotted solution are such that 1.5 < x\ + x\ < 30.817. 
Notice also that the computed solution is trapped in a neighborhood of T, outside 
7, when E is not attractive . Periodie motion persists also for smaller values of t, 
though the orbit shrinks to the origin as r —)• 0 . 

d) Unregularized integration. Here, this approach works well for relaxed values of 
the tolerances. To witness, the numerical solution computed with ode23s and 
RelTol=AbsTol=10“^, stays elose to E up to x\ + x\ ~ 2.0014. It then leaves 
E to slide on Ej". With lower values o/RelTol and AbsTol, the integrator takes 
more than half an hour in the time interval [0, 0.5]. The integrator ode45 shows a 
similar behavior. 


Example 4.2 (E looses attractivity through a non tangential potential exit point). The 
vector fields are 


/( 3 -X 3 )/ 5 \ 

(4.2) /i =1 -1/5 1 , /2 



h = 



and E is the x^-axis, with uniquely defined sliding vector field X3 = 1. E is attraetive 
through sliding along Ej*" and E^ for X3 < 3. At x^ = 3, /i is tangent to Ei and points 
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Figure 3. Example 14.11 Numerical solution computed with “Determinis¬ 
tic Euler” and constant stepsize r = 10“^. 



Eigure 4. System ()4.2p . Average trajectory computed with Random Eu¬ 
ler and stepsize r = 10 “®. 



away from S so that the point (0,0,3) is a potential non tangential exit point. When 
X 3 > 3, the vector field fi points away both from Si and S 2 and we expect the solution of 
(jl.ll) to move away from S and to enter Ri. 

a) Random Euler. This behaves qualitatively as we expect. We take 100 initial con¬ 
ditions uniformly distributed in [—r, r]^. With r = lO”'^, the average “exit point” 
is X = 2.9854 with standard deviation Xg = 0.0122. With r = 10“®, we obtain 
X = 2.9923, Xg = 0.0046. The average trajectory obtained with stepsize r = 10“® 
is plotted in Figure^ The exit point for the trajeetory is x zz 2.9944. 

b) Regularized Integration. For the regularized vector field (j2.2l) . we take a and /3 as 

in but we make different choices for the parameter e. We use two different 

parameters for a and (3, namely Ca = ep = 10“® and Ca = 10“®, ep = 10“®. The 
corresponding fast systems (|2.8I) and (12.9|) are not orbitally equivalent and, as a 
consequence of this, the behavior of the regularized solutions differs as well. We 
first consider the regularized system for Ca = £/?• In this case, {a* (xs), f* (xs)) 
is an asymptotically stable equilibrium of dlSD in (0,1)2 for X 3 < X 3 3.3853. 
For X 3 = X 3 , the eigenvalues of the Jaeobian matrix of dlSl) at the equilibrium 
cross the imaginary axis and {a*(x^),/3*(x^)) undergoes a subcritical Hopf bifurea- 
tion. In Figure\^ on the left, we plot both the stable equilibrium and the unstable 
periodic orbit for (12.8h and X 3 = 3.8. The periodic orbit separates the solutions 
that converge to the equilibrium from the solutions that reach (possibly in 
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Figure 5. Example 14.21 Left: unstable periodic orbit and equilibrium 
(a*,/3*) of (j2.8j) for xs = 3.8. Right: unstable periodic orbit and equilib¬ 
rium of (I2.10j) for X 3 = 3.447. 



finite time) the boundary o/[0,1]^. Notice that, if the initial condition of (|2.4I) 
satisfies (Hi) of Theorem \2.3l then b) in Theorem \2.3\ follows, and the regularized 
solution remains in a neighborhood of T, as long as I3*{xfi)) is attractive. 

As a consequence of this reasoning, for e sufficiently small, we expect the solution 
of (j2.4p to remain close to S up to xs = x^. For x^ > x^, all solutions inside 
( 0 , 1 )^ reach the boundary of the square, and we expect the regularized solution 
to move away from S. This is well illustrated in the left plot of Figure^ The 
equilibrium of (12.101) undergoes a subcritical Hopf bifurcation at 0:3 = xs ~ 3.4476. 
For X 3 > X 3 , all solutions leave the square [0,1]^. In the right plot of Figure\^ we 
show the unstable periodic orbit together with the stable equilibrium for (I2.10p with 
X3 = 3.447. 

With the second set of parameter values, = 10“^ and ep = 10“®, the dynam¬ 
ics of ([22]) is similar, but the bifurcation values are different. The equilibrium 
{a* {X 3 ), P* {X 3 )) of (12.9|) is stable up to X 3 = X 3 ~ 3.2296 and at X 3 = X 3 it un¬ 
dergoes a subcritical Hopf bifurcation. In the left plot of Figure 0 we plot both the 
unstable periodic orbit and the stable equilibrium of (j2.9p for ep = 10“^ 

and X3 = 3.22. We also integrate the fast system backward in time for X3 = 3.21 
from an initial condition in a neighborhood of the equilibrium. There is no periodic 
orbit in this case, since the backward solution meets the boundary o/[ 0 , 1 ]^ in finite 
time. The star in the plot marks the equilibrium (a“,0) of the fast system. This 
is the sliding vector field on 

On the left of Figure we plot the numerical solution obtained integrating the 
regularized problem (|2.2I) . X = f%{x), with ode23s and RelTol=AbsTol= 10 
for Ca = cp = 10“®. The continuous line is the first component of the solution, 
while the dotted line is the second component. They are plotted against the third 
component. The solution stays close to S up to X 3 ~ 3.39 and then it enters 
Ri. This is in agreement with the singular perturbation analysis. For the second 
set of parameters, e^ = 10 “®, ep = 10 “^, the numerical solution obtained with 
ode23s remains close to S up to X 3 ~ 3.369, and this is not in agreement with 
the singular perturbation analysis. The numerical solution obtained with odel5s 
is even less accurate and exits at X 3 ~ 3.485. We integrated the problem also 
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Figure 6 . Example 14.21 Left Plot; unstable periodic orbit and stable 
equilibrium for the fast system ( 12 .9p . = 10 " e /3 = 10 ® and xs = 3.22. 
Right plot: backward trajectory for the fast system for X 3 = 3.21. 



Figure 7. Example 14.21 Solution of the regularized system (12.2p with 
ode45. Left: = 10“®. Right: Cq, = 10“® and ep = 10“®. 




with lower values of e^, while still having €a = -^ so that the corresponding fast 
systems are all equivalent. The numerical approximations are even less reliable, 
and indeed the computed exit values increase as eg decreases. On the other hand, 
the numerical solution computed with ode45 behave more reliable leaves S to enter 
Ri for X3 ~ 3.238, as we can see from the right plot in Figure^ 


Example 4.3 (Spiral dynamics around S). In this example the vector fields are 



T, is the xs-axis, with uniquely defined Filippov sliding motion X 3 = 1, and there is spiral 
like dynamics around S. For X 3 < 1, T, is attractive in finite time while, for X 3 > 1, Tj is 
not locally attractive. 

a) Random Euler. Approximations computed with “Random Euler” move away from 
S for X3 > 1. In FigurelEwe show the average trajectory computed with r = 10“^ 
and initial condition (10“®, 10“®, 0.5). For a statistical estimate of the exit point, 
we consider an ensamble of 100 initial conditions with the first two components 
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Figure 8 . Example 14.81 Average trajectory obtained with Random Euler 
and r = 10 “®. 


x10' 



Eigure 9. Example I4.3I Left: a solution of (|4.3p obtained with Random 
Euler and r = 10~®. Right: norm of the vector (hi,/i 2 ) in function of X 3 
for this trajecytory. 




uniformly distributed in [—r, r]^. For r = 10“^, the mean value of X 3 at the exit 
point is X 3 ~ 0.94777, with standard deviation ~ 0.04297. For r = 10“®, the 
mean value of X 3 at the exit point is X 3 ~ 0.97910 and the standard deviation is 
~ 0.00529. In Figurel^ on the left, we plot a trajectory of (HD obtained with 
Random Euler and stepsize r = 10“®. On the right of Figure\^ we plot the 2-norm 
of the vector (/ii,/i 2 ) (for us, this is the vector {xi,X 2 )) in function of X 3 . The 
“X” in the plot marks the estimated exit point. 
b) Regularized Integration. Next, con.sider integrating the regularized vector field (| 2 . 2 |) 
with a and ft as in (12.3p and two different sets of parameters. We first consider 
Cq = e /3 = 10“^. For these parameters values, the equilibrium {a*, (I*) of the 
fast system (12. 8 |) is stable up to X 3 ~ 1.41798, then the eigenvalues of the Jacobian 
matrix at the equilibrium cross the imaginary axis. We compute the solution of the 
regularized system with initial condition (10“^, 10““^, 0) in the time interval [0,2] 
(so thatO <X 3 < 2 ) with ode23s, odelBs and ode45 raif/i RelTol=AbsTol= 10“^^. 
The computed solutions show different behavior: the numerical solution computed 
with ode23s and odel5s remain in a small neighborhood of S for the whole time 
interval [0,2], while the solution computed with ode45 moves away from E for 
X 3 ~ 1.5 as it is seen from the left plot of Figure fllK The average stepsize used 
by the stiff integrators is 0(10“^), while the one used by the explicit integrator is 
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Figure 10. Example 14.31 Regularized integration with ode23s. Left; 
= 10“^. Right: ea = 10“^ and ep = 10“^. 




Figure 11. Example 14.31 Unregularized integration with odel5s and 
RelTol=AbsTol=10"®. 



0(10 ^). We then consider a second set of parameters = 10 ^ and = 10 
The equilibrium {a*, 13*) of the corresponding fast system (EJ]) is unstable when 
X 3 > 0. On the right of Fiaure\l(]\ we plot the approximation computed with ode23s 
(the approximation computed with ode45 behaves in a similar way). 
c) Unregularized integration. The numerical solution obtained with the stiff Matlab 
integrator ode23s, with RelTol= AbsTol = 10“^, stays close to T, up to x^ = 1 
and then it leaves T, to enter one of the Rj’s. For lower values of RelTol, the 
integrator takes more than half an hour in the time interval [0,1]. The numerical 
solution obtained with odel5s and RelTol=AbsTol=10“® is very inaccurate as 
it is evident from the plot in Figure \4-3[ Lower tolerances do not produce better 
results. 
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Figure 12. ExamDle l4.4l plot in the plane (xs, X 4 ) of 1000 approximations 
obtained with Random Enler and r = 10“^, and initial condition 
( 10 - 2 , 10 - 2 , 3 , 0 . 9 ). 



Example 4.4 (Filippov’s vector field on S is ambiguous). We consider the following 
vector fields 


(4.4) 


( ^ ^ 

f = ^ 

^ -X3 + ix4 ’ 

V X4 / 

/-(X3 - 3)2 - (X4 - 3)2 + 5\ 

1 

— X3 + 28X4 

V X4 / 


1 

2 1 

-X3 + 23^4 

V X4 / 

-1 

-X3 + 4X4 ’ 

V X4 / 


and S is the (X 3 , X 4 ) plane. The circle 7 = {x E (X 3 — 3)2 + (x 4 — 3)2 = 4} (see Figure 
CD, divides E in two regions. Outside 7 , E is attractive upon sliding along E^ and E^. 
For all points on 7 , the vector field is tangent to E and, inside 7 , it points away from 
E, so that E is not attractive from inside 7 . We would expect the solution of (dl) to 
leave E once it reaches 7 . Note that there is a family of Filippov sliding vector fields on 
E, namely: xi = X 2 = 0, X 3 = —X 3 + (16 — 13A)x4, X 4 = X 4 , with 0 < A < |. 


(a) Random Euler. The ambiguity of a Filippov sliding vector field is clearly reflected in 
the numerical solutions computed with Random Euler. In Figure [221 we plot the X 3 
and X 4 eomponents of 1000 trajectories computed with r = IO- 2 , and same initial 
condition (10-2,10-2,3,0.9). The dotted cirele in the plot is the eurve 7 . The 
two bold darker lines are two sample trajectories. The shaded region is obtained by 
plotting all 1000 trajectories. The plot suggests that the ehoiee of random stepsizes 
covers the region obtained by choosing one of the possible vector fields in Filippov’s 
differential inclusion. For completeness, in Figure [23 we plot in function of time 
the first and second component of the average trajectory computed with Random 
Euler with T = 10-^. Att'zs. 0.1067, the third and fourth eomponents of the average 
trajectory are on 7 and indeed the plotted solution leaves E and starts sliding on 
E^ at t zz 0 . 12 . 
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Figure 13. Example 14.41 First and second component of the numerical 
solution computed with Random Euler with r = 10“®. 



b) Regularized Integration. Here, the computed approximations of the regularized sys¬ 
tem move away from S once they reach 7 . We note that j is a curve of saddle- 
node bifurcation values for (USD for any parameter value and/or ep. On 7 , 
{ a *,/*) = (1,0.5) is a double root of (j2.7p . while inside 7 there is no solution of 

(EZl). 


5. Conclusions 

In this work we have been interested in studying the behavior of solutions of piecewise 
smooth systems in the neighborhood of a co-dimension 2 discontinuity surface S, intersec¬ 
tion of two co-dimension 1 discontinuity surfaces. It has long been accepted that if solution 
trajectories cannot leave S (S is attractive), some form of sliding motion on E should be 
taking place. Precisely which sliding motion has been the subject of much investigation, 
but it has not been our concern in this paper. Our chief interest in this work has been 
trying to understand what should happen when S loses attractivity (at generic first order 
exit points). To our knowledge, this type of study had not been carried out before. 

We took the point of view that the piecewise smooth system (HI) was the only informa¬ 
tion at our disposal, and treated this model with its own mathematical dignity. Naturally, 
if (HI]) arose as a simplified model for some other known differential system, then this 
original system should ultimately guide the search for appropriate dynamics near S, and it 
may well be that the dynamics of this “true” system are not matched by those of (jl.ip . If 
this is the case, we should legitimately question the validity of the model ()l.ll) in the first 
place. On the other hand, in the absence of knowledge of an underlying “true” system, 
when (jl.ip is the only datum we have, then we believe that we should try to modify this 
model so that the dynamics of the modified system match those of (jl.ip . 

To obtain information on the dynamics of we proposed use of a simple Euler 

method with random steps, uniformly chosen with respect to a reference, small, stepsize. 
Our study unambiguously show that: (i) when E is attractive, solution trajectories remain 
near E (thereby validating an idealized sliding motion on E); (ii) when E loses attractivity, 
solution trajectories leave a neighborhood of E. 
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Several other possibilities have also been considered in this work: regularization tech¬ 
niques, plain and simple Euler method with fixed stepsize, and direct numerical integration 
of (jl.ip with sophisticated off-the-shelf solvers for differential equations. None of these 
options satisfactorily resolved the dynamics of (ll.lh . and often produced misleading be¬ 
havior. Ultimately, this occurred because each of these choices either superimposed its 
own dynamics on those of (II.Ih (as Euler method and regularization techniques do, further 
producing different behaviors depending on how the regularization is made), or just failed 
to produce reliable answers in too many cases (this was the case with directly solving (II.Ij) 
with existing software, where the outcome dramatically dependend on the solver used, or 
on the tolerances values, or both). 

Unfortunately, our conclusions are not fully satisfactory either. Our analysis tells us 
that is the dynamics of O) around S that must be used to tell us what should happen in 
a neighborhood of S, but we know of no general foolproof mean to regularize the system 
so that the regularized trajectory will be following the dynamics of (jl.ip . Perhaps, and 
-again- as long as the model is appropriate, the most reliable and practically efficient 
way to proceed is to accept some form of idealized sliding motion on S as long as S is 
attractive, while also demanding that a sliding trajectory leaves S when the latter loses its 
attractivity. The construction of appropriate sliding vector fields fulfilling these requests 
remains an outstanding and challenging task. 
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